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Abstract. We discuss a new kind of astrophysical transport problem: the 
coherent evolution of neutrino flavor in core collapse supcrnovac. Solution 
of this problem requires a numerical approach which can simulate accurately 
the quantum mechanical coupling of intersecting neutrino trajectories and the 
associated nonlincarity which characterizes neutrino flavor conversion. We 
describe here the two codes developed to attack this problem. We also describe 
the surprising phenomena revealed by these numerical calculations. Chief 
among these is that the nonlincaritics in the problem can engineer neutrino 
flavor transformation which is dramatically different than in standard Mikheyev- 
Smirnov-Wolfenstein treatments. This happens even though the neutrino mass- 
squared differences are measured to be small, and even when neutrino self-coupling 
is sub-dominant. Our numerical work has revealed potential signatures which, if 
detected in the neutrino burst from a Galactic core collapse event, could reveal 
heretofore unmeasurable properties of the neutrinos, such as the mass hierarchy 
and vacuum mixing angle #13. 



PACS numbers: 02.70.-c, 14.60.Pq, 97.60. Bw 
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1. Introduction 

The ghost-like neutrinos (v e , P e , v^, u^, v T and D T ) are chargeless, spin-| particles 
notorious for their weak interactions. The revelation by recent observations and 
experiments of non-zero neutrino rest masses and vacuum neutrino flavor mixing 
has forced the astrophysics community to confront a vexing problem: the nonlinear 
evolution of neutrino flavor in astrophysical environments where there are appreciable 
neutrino fluxes. Solution of this problem could be important because many 
environments associated with compact objects and the very early universe are 
energetically and in other ways dominated by neutrinos and their interactions. Perhaps 
the most impressive example of this is in core collapse supernovae, where the collapse 
of a highly-evolved core supported by relativistically-degenerate electrons to a cold 
neutron star releases ~ 10% of the core's rest mass in neutrinos of all kinds (see, e.g., 
[1] for an introductory review and [2] for a discussion of general supernova neutrino 
physics issue). 

Neutrinos in this environment are thought to play a pivotal role in nearly every 
aspect of supernova physics, from the explosion mechanism itself, at fairly early times 
after the collapse of the core, to setting the conditions for the synthesis of heavy 
elements at later times. By virtue of their weak interactions and tenuous coupling to 
ordinary matter, neutrinos can transport energy, entropy, and lepton number through 
very dense matter that other particles might not be able to penetrate. And neutrinos 
can more than make up for their feeble individual interactions with huge numbers. 

The astrophysics community has long recognized the important role that 
neutrinos play in core collapse supernova explosions. From the earliest work of Colgate, 
Wilson, and others [3, 4, 5, 6] to the most recent and most sophisticated numerical 
simulations (e.g., [7, 8, 9, 10, 11, 12]), the key problems have been to characterize 
the transport of neutrinos in and above the neutron star core and to calculate the 
energy spectra and fluxes of each neutrino species. Solving these problems has required 
not only front-line coupled multi-dimensional hydrodynamics and Boltzmann neutrino 
transport computations but also understanding physics which is not readily accessible 
in the laboratory, such as the equation of state of hot, neutron-rich nuclear matter. 

However, this game changes somewhat in light of the experimental fact of neutrino 
flavor transformation. It has been shown that v e and D ei neutrinos and antineutrinos 
with electron flavor (i.e., neutrinos associated with electrons and positrons in weak 
interactions), can be transformed into neutrinos and antineutrinos with other flavors, 
i/„, v T , and v r (see, e.g., [: ] for a review). Now both goals of calculating 
neutrino transport and finding the emergent neutrino energy spectra and fluxes may 
not be attainable without an adequate treatment of medium-affected neutrino flavor 
transformation. As we shall see below, in many circumstances adding on a flavor 
evolution solver would greatly complicate neutrino transport calculations. We can 
identify two broad regimes in this problem: (1) the high density partially coherent or 
non-coherent regime, where the neutrino transport mean free paths are comparable to 
or smaller than neutrino flavor oscillation lengths or Mikkheyev-Smirnov-Wolfenstein 
(MSW) [14, 15, 16] resonance widths; (2) the coherent regime, where neutrinos are 
essentially free streaming on scales relevant for neutrino flavor transformation. The 
first regime can be very difficult to treat and requires a quantum kinetic approach 
(e.g., [17, 18, 19, 20, 21, 22]). However, it is generally thought that active-active 
neutrino flavor transformation is negligible in this regime. This prejudice is motivated 
by the fact that the relevant active-active channel neutrino mass-squared differences 
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Figure 1. Illustration of the quantum coupling problem for neutrinos on 
different, but intersecting trajectories. The world lines for three neutrino beams, 
Vk, v p and v q , intersect as shown at points O, P and Q. Neutrino-neutrino 
forward scattering at these intersection points will quantum mechanically couple 
the flavor evolution histories of these neutrinos. 



are small, the matter density is high in this regime [23, 17], and collective neutrino 
flavor transformation is suppressed in very dense neutrino gases [24]. Our numerical 
approach concentrates on the second regime. 

The coherent neutrino flavor evolution problem differs in two essential aspects 
from conventional transport problems. First, of course, is quantum mechanical 
coherence itself. Second is neutrino self-coupling. 

Coherence dictates that we follow the development in time and/or space of the 
complex amplitudes which describe the flavor states of neutrinos. This must be done 
with a Schrodinger-like equation. On a given neutrino trajectory, or world line, with 
Affinc parameter t, this is 



Here the operator H is a Hermitian Hamiltonian operator that contains the vacuum 
and in-medium generators of neutrino flavor evolution, while \ip (t)) is the ket that 
describes the flavor state of a neutrino at point t on this given trajectory. Neutrinos 
in supernovae for the most part will be ultra-relativistic, so it may seem odd that 
a "non-relativistic Schrodinger equation" can give an adequate account of flavor 
evolution. It works for two reasons [25]: (1) we consider only neutrino forward 
scattering interactions and neglect inelastic scattering; and (2) in the ultra-relativistic 
limit the information in the four-component Dirac spinor describing the neutrino's 
flavor state is mostly contained in two components, the two "large" components. 

A serious complication in following neutrino flavor transformation in supernovae 
stems from neutrino-neutrino forward scattering contributions to H. These render 
the problem nonlinear, since the interactions which dictate flavor transformation 
amplitudes are themselves dependent on the neutrino flavor states. This seemingly 
innocuous statement masks a kind of geometric coupling which, to our knowledge, 
is unique in astrophysical transport problems. Neutrino-neutrino scattering could 
quantum mechanically couple the flavor evolution on all intersecting trajectories. This 
is illustrated in Figure 1. Following the flavor evolution of neutrino v p would entail 
knowing the flavor states of all neutrinos on which it forward-scattered. That might 
include, for example, neutrinos Vk and v q propagating on world lines which intersect 
Vp's trajectory at points Q and O, respectively. Note, however, that z^'s and v q 's flavor 
histories are not independent, since they may have undergone a forward scattering 
event at point P. Obviously, with a large number of neutrinos this quickly can become 
a daunting problem! 

To simplify the problem we adopt the spherically symmetric "neutrino bulb" 
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model (hereafter the "bulb model") [26]. In this model neutrinos are emitted half- 
isotropically from a spherical "neutrino sphere" (essentially the hot neutron star 
surface) with energy spectra that are initially of black body form. We adopt a smooth 
matter density profile p which depends only on radius r, the distance from the center 
of the neutron star. At late times, close to the neutron star where matter tends to be 
isothermal, this matter profile falls exponentially. Further out where the envelope may 
be closer to constant entropy in character, the density profile is p ~ r -3 . Because of 
the assumed spherical symmetry, neutrino trajectories with the same emission angle 
$o become equivalent, and neutrino modes emitted in the same flavor and energy state 
and with the same emission angle have identical flavor evolution history [26]. Here 
the emission angle $o is the angle between the neutrino propagation direction and the 
vector normal to the neutrino sphere at the emission point. 

In the bulb model, all neutrino evolution above the neutrino sphere is assumed 
to be coherent and essentially free-streaming in character. This will be a good 
approximation for the conditions expected to obtain a few seconds after core bounce, 
and may be useful for earlier times in some cases. We bin neutrino modes emitted in 
each flavor by both the emission angle and the energy of the neutrino mode. Because 
the flavor histories of all neutrino bins are coupled, we solve self-consistently for the 
flavor evolution of all neutrino modes simultaneously. The flavor evolution of each 
neutrino mode is described by an equation with the form of Equation (1). Each 
neutrino trajectory has an individual Afline parameter which can be expressed in 
terms of radius r. 

There are a variety of physical processes which can affect neutrino flavor evolution, 
either by complicating the coherent calculation procedure outlined above, or by 
introducing decoherent effects. As an example of the former, the density profile in 
the core collapse supernova environment can be quite dependent on time, with shocks 
and multi-dimesional effects of paramount importance in some epochs (e.g., [27, 28]). 
In the latter case, matter density fluctuations associated with shocks, sound waves, and 
turbulence can lead to neutrino flavor depolarization and decoherence [29, 30, 31, 32]. 

In our numerical work we have chosen to concentrate on the simplest environments 
in order to gain insight into the effects of coupled quantum mechanical flavor evolution. 
Our rationale is that this physics will be present always in the environments we treat. 
The nonlinear collective phenomena revealed by our calculations can be modified in 
real supernovae by the kinds of effects discussed above. However, our calculations lead 
us to conclude that most of the qualitative and even some of the quantitative results 
for our simple models are likely to survive. Even our simplified supernova models 
present novel numerical challenges which will have to be faced in any comprehensive 
treatment of the subject. 

2. Solution method 

Even with the simplified model discussed in Section 1, obtaining a convincing solution 
to the full, coupled-trajectory supernova neutrino flavor evolution problem, which we 
term the "multi-angle" problem, is still a challenging task. So difficult, in fact, that 
there had been no attempt to solve it prior to our effort. This meant that previous 
work provided no guidance on numerical strategies for solving multi-angle neutrino 
flavor evolution. The only clues to general coherent neutrino flavor phenomena in 
supernovae came from so-called "single-angle" calculations, where trajectories were 
coupled, but the flavor evolution on every trajectory was taken to be the same as 
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that on a radially-directed trajectory. In Section 2.1 we outline the plan we used to 
attack this problem. Partly because of the lack of previous work, we have developed 
two numerical codes more or less independently at UCSD and LANL. As an aid to 
debugging and as a device for teasing apart physics issues from numerical issues, the 
results from the two codes were compared for key test problems. This procedure 
proved efficacious, in part because sometimes numerical issues and problems tended 
to be different in each code. In Sections 2.2-2.4 we highlight the key implementations 
of FLAT, the computer code developed at UCSD. In Section 2.5 we discuss some 
features of BULB, the code developed at LANL. 

2.1. Project plan 

Our ultimate goal is to solve for the flavor histories of neutrinos and antineutrinos (i.e. 
v e , Vn, v T , v e , Va and D T ) emitted in the bulb model with realistic energy spectra and 
traveling on all physically realizable trajectories. This is what we term the "three- 
flavor multi-angle" problem. 

In the literature, the neutrino oscillation problem with three flavors is 
usually approximated as two two-flavor oscillation problems that occur at different 
time/distance scales. The rationale for this is that the neutrino mass-squared 
differences associated with these two-flavor oscillation problems are different by more 
than an order of magnitude. Previous work on medium-affected neutrino oscillations 
exclusively adopted this two-flavor approximation [17, 33, 34, 35, 36, 37, 38, 39]. 

These studies also adopted what we termed above the "single-angle" 
approximation. In this approximation flavor evolution of neutrinos with energy E 
propagating along different neutrino trajectories is assumed to be the same as that 
of neutrinos with the same energy E propagating along a representative trajectory, 
usually taken to be the radial direction. The validity of both the two-flavor and single- 
angle approximations for supernova neutrino flavor evolution, however, remains to be 
shown. 

Because of the big gap between the ultimate goal we would like to achieve and 
the state of art of neutrino oscillation studies, it is clearly very risky to attempt a 
solution to the full problem in one single step. Not only will there be little theoretical 
guidance in such a monolithic approach, but also finding physical interpretations of 
any novel numerical results could prove to be problematic. Instead, we have sought 
to achieve our goal in four stages. 

In the first stage we solve for neutrino flavor transformation using the two-flavor, 
single-angle approximation. This stage overlaps with studies which already exist in 
the literature. This procedure should offer a good check on our numerical codes as 
well as on existing theories. 

In the second stage we will still assume two neutrino flavors, but will solve for 
neutrino and antineutrino flavor evolution using a full multi-angle implementation of 
the bulb model. The difference between this stage and the first stage should give us 
insights on how multi-angle effects may enter into neutrino flavor transformation and 
how good the single-angle approximation might be. 

In the third stage we will implement full three-neutrino flavor evolution but with 
the single-angle approximation. Comparison of the results of this stage and those of 
the first stage should give us important clues on whether the two-flavor approximation 
is valid in collective neutrino oscillations. 

In the last stage we will solve the full problem, i.e., neutrino oscillations with 
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Figure 2. Illustration showing the structure of FLAT. This code is a collection 
of modules. These modules are represented here as various building blocks in 
the left part of the figure. An executable binary can be compiled by choosing 
appropriate modules for a particular problem. This is illustrated in the right 
part of the figure. Building blocks with the same geometric shape and the same 
top/bottom layer(s), but with different colors, represent modules that have the 
same interface functions but perhaps employ different physical approximations 
or numerical algorithms. Changing a particular physical setting or a numerical 
algorithm can be achieved by swapping out the appropriate module(s) in the code 
with other module(s) with the same interface functions. 

three neutrino flavors and a full multi-angle implementation. Presumably, the results 
of the computation in this stage could be estimated at least crudely from the theories 
which employ the two-flavor, single-angle model and the results of the second and 
third stages calculations. Comparison between the numerical results in the full three- 
neutrino calculations and the theoretical estimates again can be used to check both 
the numerical codes and the theories. 

2.2. FLAT's highly modular design 

As explained above, we plan to achieve our ultimate goal progressively, essentially by 
solving three intermediate problems before attacking the full problem. These problems 
differ only in the physical approximations that they employ. Written in C++, FLAT 
takes advantage of the Object-Oriented Programming paradigm and provides a unique 
and uniform solution to all these problems by offering a hierarchical set of problem- 
specific modules. Changing a particular physical setting or a numerical algorithm can 
be done by swapping out appropriate module(s). The structure of FLAT is illustrated 
in Figure 2. 

At the lowest level of the hierarchy are the NeuBin modules. Each of the NeuBin 
modules implements two C++ classes, NeuPot and NeuBin, which correspond to the 
Hamiltonian H and the flavor state \ip) in Equation (1), respectively. We note that 
the Schrodinger Equation (1), as well as the numerical algorithm that FLAT uses 
to solve this equation, do not depend on the specifics of the representations of H or 
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Therefore, the higher- level modules in FLAT are independent of the internal 
implementation in the NeuBin modules and can operate through a uniform interface 
on any NeuBin module. The NeuBin modules are represented as pentagonal blocks in 
Figure 2. The NeuBin module interfaces are represented by the common-color layers 
on top of the pentagonal blocks. 

Two examples of the interface functions of a NeuBin class are NeuBin: : evolve () 
and NeuBin: :diff (). Function NeuBin: : evolve () evolves the corresponding flavor 
state \ip) residing at Afhne parameter value ton a given trajectory one step further for 
given Hamiltonian H and step size At (see Section 2.3). Function NeuBin: :diff () 
computes the "difference" between two flavor states. This difference is used to estimate 
numerical errors. 

There are two NeuBin modules that have been used in production runs: 
NeuBin_F2C and NeuBin_F3C. These modules define the flavor state and Hamiltonian 
for two-flavor and three-flavor problems, respectively. Once the program has been 
proven to work well using the two-flavor approximation, we can use the same 
program to solve three-flavor problems by swapping out the NeuBin_F2C module with 
NeuBin_F3C. 

Above the NeuBin modules are the NeuBeam modules, each of which defines 
a NeuBeam class. A NeuBeam object contains an array of NeuBin objects. Each 
element of this array corresponds to a flavor state \ip{E)) for a particular energy 
bin [E — hAE,E + ^AE\. All NeuBeam modules share the same interface through 
which higher-level modules can operate NeuBeam objects. 

NeuBeam: :sum() is one of the important interface functions for the NeuBeam 
modules. This function does a numerical integration over all energy bins. This 
summation represents an important step in producing the neutrino self-coupling 
Hamiltonian H vv (see Section 2.5). The algorithm used in this numerical integration 
is specific to each module and, again, is hidden from higher-level modules. 

The NeuBeam module used in production runs is called NeuBeam_SS. In this module 
NeuBeam: : sum() employs a simple midpoint rule for numerical integration. If a more 
sophisticated algorithm is required, a new NeuBeam module can be created. The 
NeuBeam_SS module can be swapped out with this new module without modification 
to the remaining components of FLAT. The corresponding upgrade of the code is done 
automatically. 

On top of the NeuBeam modules are the NBGroup modules. There are three 
kinds of NBGroup modules: NBGroup_EM, NBGroup_I0, and NBGroup_PM. These modules 
specify the physical environment, the input/output (I/O) methods, and the numerical 
algorithm to solve differential equations, respectively. These three types of modules 
share the same core data, i.e., the current flavor states of all neutrinos, and, at the 
same time, perform nearly independent tasks. This special relation among NBGroup 
modules is realized through the inheritance of classes and virtual functions. An 
NBGroup class is defined first. This class contains all the data that are used by at 
least two kinds of NBGroup modules. This class also declares a variety of virtual 
functions such as NBGroup: :set_bg(), NBGroup: :check() and NBGroup: : evolve (). 
Each of these virtual functions is defined in a child class of NBGroup, but may be called 
by other child classes. 

Each of the NBGroup_EM modules defines an NBGroup_EM class as a child class 
of NBGroup. A key interface function of these modules is NBGroup_EM: :set_bg() 
which, by definition of a virtual function, is the "real function" to be executed 
when virtual function NBGroup : :set_bg() is called. NBGroup_EM: :set_bg() calculates 
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-ffmatt, the ordinary matter background contribution to the total Hamiltonian 
H. NBGroup_EM: : set_bg() also calculates -ff„„, the background neutrino medium 
contribution to H. 

There are two NBGroup_EM modules used in production runs: NBGroup _PNS1D and 
NBGroup _PNS2D. In NBGroup _PNS ID the single-angle approximation is implemented, 
while in NBGroup_PNS2D the multi-angle implementation is employed. 

Each of the NBGroup_I0 modules defines an NBGroup_I0 class as a child class 
of NBGroup. An NBGroup_I0 class defines NBGroup_I0: :check(), which is the real 
function corresponding to NBGroup :: check () . NBGroup_I0 : : checkO performs a 
variety of critical checks. For example, if necessary, it can save a snapshot of the 
current flavor states of all the neutrino modes, e.g., at specified time intervals. Also, 
NBGroup_ID: :check() can set a stop flag if the specified time allocation is used up. 

The NBGroup_I0 modules used in production runs are NBGroup_NC and 
NBGroup_NC_SMP. Module NBGroup_NC reads and writes data in the NetCDF format 
[40], a cross-platform binary data format. Module NBGroup _NC_SMP is similar to 
NBGroupJJC, except that it handles I/O in a different POSIX thread (Pthread) [41]. 
As a result, NBGroup_NC_SMP is able to complete NBGroup_ID: : check () function calls 
promptly. This helps to augment the overall code performance, because the CPU 
can continue to do computation during the time in which it would otherwise be idle 
waiting for slow I/O data transfers to complete. 

Each of the NBGroup_PM modules defines an NBGroup_PM class as a 
child class of NBGroup. An NBGroup_PM class defines a member function 
NBGroup _PM :: evolve () , the real function corresponding to NBGroup :: evolve () . 
Once called, NBGroup_PM: : evolve () can continuously evolve neutrino flavor states 
over radial distance until stop flags are encountered. 

The NBGroup_PM modules used in production runs are NBGroup_MM and 
NBGroup _MM_SMP. In the NBGroup _MM module, NBGroup _PM: : evolve () employs 
an adaptive step-size algorithm to solve Equation (1) for each neutrino 
trajectory (see Section 2.4). During each sub-step of the algorithm's execution, 
NBGroup_PM: : evolve () calls NBGroup :: set_bg() to set H mat t + H uv . At the 
end of each whole step in the algorithm's execution, NBGroup _PM :: evolve () calls 
NBGroup: : checkO to check and see if any I/O needs to be done and if the program 
needs to be terminated. NBGroup_MM_SMP implements a multi-threaded version of the 
algorithm used in NBGroup _MM. FLAT automatically uses the multi- CPU/core feature 
of a modern workstation when the NBGroup _MM module is replaced by NBGroup_MM_SMP. 

At the highest level, a Driver module defines a class designated as Neutrinos. 
This class inherits the NBGroup_EM, NBGroup _PM and NBGroup_I0 classes. In the Driver 
module, Neutrinos: : evolve () (which is actually NBGroup_PM: : evolve ()) is called 
to solve the neutrino flavor transformation problem. Driver also handles command 
line options. These might include, for example, setting the time allocation and/or the 
final radius. 

FLAT'S highly modular framework endows it with great extensibility. This 
aids program development as well as helps with the exploration of neutrino flavor 
conversion physics. With this modular architecture it is straightforward to pick 
the simplest numerical algorithm to attack a given problem. Additional or different 
modules can be added easily if more sophisticated algorithms are necessary, or if the 
intention is to treat different problems. When following the project plan outlined 
in Section 2.1, the modular framework of FLAT enables us to minimize redundancy 
in code writing and maximize usage of existing code. It also maintains consistency 
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when we progress from one stage to the next and, therefore, isolates the differences 
in various physics approximations when we compare the results obtained in different 
stages. 



2.3. NeuBin modules in FLAT 

Currently FLAT can utilize two different NeuBin modules: NeuBin_F2C and 
NeuBin_F3C. These implement the Hamiltonian H and the neutrino flavor state \ip) 
for the 2-flavor and 3-flavor mixing schemes, respectively. 

The NeuBin_F2C module defines two classes, NeuBin and NeuPot. A NeuBin 
object contains two complex variables. These variables correspond to the two complex 
components of the flavor wavefunction 

« ,s («:)• < 2 > 

where a Va = (v a \ip) is the amplitude for a neutrino in state to be in flavor state \v a ) , 
and where the flavor index a = e, x, with x representing either the /iorr neutrino or a 
linear combination of these. Class NeuPot, however, does not render the corresponding 
flavor-basis representation of the Hamiltonian H as the familiar Hermitian matrix in 
the Schrodinger equation. In the 2-flavor mixing case, for example, the Hamiltonian 
would be a 2 x 2 complex array 

/ {v e \H\v e ) {v e \H\v x ) \ 

\ ( Vx \H\v e ) {v x \H\y x ) J' 1 ' 

Instead, NeuPot stores the Hamiltonian in a three-component real vector h which is 
defined by 

H = h + h-(T. (4) 
Here the three Cartesian components of the vector a corresponds to the Pauli matrices: 

ffl = ( ? o ) ' 0-2 = ( i o 1 ) ' and CT3 = ( -°i ) • (5) 

The trace of matrix H [h in Equation (4)] is of no physical significance and is 
ignorable. This is because ho generates only an overall common phase in the neutrino 
flavor amplitudes and, therefore, is not relevant for neutrino oscillations or neutrino 
flavor transformation. 

Equation (1) has an approximate solution for a small parameter step At, 

%j)(t + At) ~ exp(-iHAt)ip(t) (6) 

cos(hAt) — i/13 sin(hAt) —(ihi + /12) sin(hAt) 
— (ihi — hq) sin(ft,Ai) cos(hAt) + i/13 sin(hAt) 

where 



and 



m, (7) 



h = \h\ = J hf + h 2 2 + h% (8) 



l, . \ (9) 



With a specified Hamiltonian H (a NeuPot object), Equation (7) is employed in 
module NeuBin_F2C to evolve a NeuBin object by a small step At. The procedure 



Simulating nonlinear neutrino flavor evolution 



10 



implied by Equation (6) preserves the unitarity of ip and becomes exact when H is 
time independent. 

For a particular momentum mode A, the neutrino (antineutrino) density operator 
px (px) provides the necessary information to construct a neutrino-neutrino forward 
scattering background potential H vv for a given test neutrino. The density operator 
for neutrino mode A is, for example, p\ = \tp\){ip\\. Potential H vu is essentially a 
summation of densities operators over all neutrino and antineutrino modes, weighted 
by the directional dependence of the weak current-current Hamiltonian (see, e.g., [26]). 
The contents of the density matrices are stored in NeuPot objects. In the NeuBin_F2C 
module, this is done by, for example, writing density matrix p as 

P=\+S-<T, (10) 

where 

(Re(a* c a I/x ) \ 
Im(aZ.a„.) . (11) 

KkJ 2 -kJ 2 ) / 

The approach adopted in FLAT for handling the full 3x3 neutrino flavor evolution 
problem is essentially similar to the method outlined above for the 2x2 case. It is, of 
course, necessarily more complicated. Like the NeuBin_F2C module, the NeuBin_F3C 
module defines a NeuBin class implementing neutrino flavor wavefunction 




i' = <v , (12) 



in obvious analogy to Equation (2). The NeuPot class in the NeuBin_F3C module 
defines 9 real variables: «n, U22, "33, «12, U23, U31, W\2, u>23, and W31. Variables Uij 
and Wij are defined in 

Hij =Uij +\Wij, (13) 

where are the elements of the 3x3 flavor-basis Hamiltonian matrix 

/ {v e \H\v B ) {v e \H\vp) {v e \H\v T ) \ 
H= \ {u^\H\u e ) ( Vv \H\vp) ( Vll \H\i> r ) ■ (14) 
V (u T \H\u e ) {v T \H\v^) {v r \H\v T ) J 

The 3x3 module NeuBin_F3C also uses Equation (6) to evolve a neutrino flavor 
wavefunction over a short step At. However, it is difficult to express exp(— iHAt) for 
the 3x3 case in an explicit form like Equation (7). Instead, the module NeuBin_F3C 
solves the eigenvalue equations 

J2 H io V jk = CkV ik , i,j,k = 1,2,3, (15) 



and uses 



-id At 



cxp(-iffAt) = V e -K2At y t ( 16 ) 



g-i&At 



Here Vik represents the zth component of the fcth eigenvector of the Hamiltonian H, 
with corresponding eigenvalue £fc- 
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Module NeuBin_F3C employs an analytical method to solve Equation (15). This 
method is based on the algorithm developed in [42] and takes into account the fact 
that the trace term in H is not relevant for neutrino flavor evolution. Assuming that 



un + U 2 2 + u 33 = 
and, if necessary, shifting uu to make this sum 0, we define 

-3[(U22U33 - - («l 2 + W 2 12 ) - (ti| 3 + w| 3 ) - + U&) 



w ll) + U3, 3 {u\ 2 + w\ 2 ) - U11U22U33 



p 

q= -^[W11(W23+ U '23) + U 22('U31 

+ 2(-Mi 2 M 2 3'"31 + W12W23W3I + W12U23W31 + W12W23U31)], 

(j) — arccos(q). 

The eigenvalues can be written as 

2 Vp , 



(17) 

(18) 

(19) 
(20) 

(21a) 

(216) 
(21c) 

From Equation (15) it can be seen that the eigenvector Vj (with corresponding 
eigenvalue d) is orthogonal to the complex conjugates of the row vectors in matrix 
H — Qil. Assuming that X^ and X^ are two linearly independent row vectors in 
matrix H — QI, we can write 

* (1) ** (2) (22) 



3 

2VP 



cos 



2tt 

T 



C1-C2 



Vi = 



XW x X( 2 )| 



Because the eigenvectors of H are orthogonal to each other, it is necessary to solve for 
only two eigenvectors, say V± and V2. The third eigenvector, V3, can be determined 
from 



V 3 = V{ x V 2 *. 



(23) 



2-4- NBGroup_PM modules and the parallelization of FLAT 

Module NBGroup_MM is a single-thread version of an NBGroup_PM module. This single 
thread implementation uses an algorithm, similar to the modified midpoint method 
[43], to determine \ip\(r + Ar)) from |i/ , i°' l ( r )) according to Eq. (1). (Note that the 
relation between the radial coordinate r and the Affine parameter value t for each 
mode A depends on the trajectory direction of this neutrino mode.) The algorithm 
consists of three parts. In the first part, an approximation to \ip\(r + Ar)) is obtained 
using radial coordinate step size Ar: 

Step 1: Call NBGroup : : set_bg() to set the background Hamiltonian H^($q) = 
H ma ,tt{n e (r)) + H ul ,(ip^ \r,-&o) for each neutrino trajectory i?o- This function 
utilizes flavor states |V'i^( r )) f° r a ^ neutrino modes A' in the problem. 

Step 2: Find \ip^\r + Ar)) := U\(H^ g , r, Ar)\ip^\r)) for each neutrino mode A. 

Here U\(H^,r,Ar) is the operator evolving |^) through a small step. This 
"parameter evolution operator" is implemented as NeuBin: : evolve (). The 
step in mode A's Affine parameter f# po corresponding to Ar is Ai$ (r, Ar) = 
U (r + Ar)-U (r). 
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Step 3: Compute H^^o) = H m&tt (n e {r + Ar)) + + Ar,z9 ) for each 

neutrino trajectory $ . 

Step 4: Find \ip^\r + Ar)) := U\(H^} , r, Ar)|-0^ O ' l (r)) for each neutrino mode A. 

Step 5: Find (V'i ( r + ^ r )) as an average of + Ar)) an d IV^C 7 " + Ar)) f° r 

each neutrino mode A. 

The second part of the algorithm essentially repeats the first part, but with the 
step size reduced to ^Ar: 

Step 6: Find \^\r + §Ar)) := U(H^,r, |Ar)|^i 0) (r)). 

Step 7: Compute H$(#o) = H matt (n e (r + §Ar)) + ff^V + iAr,i? ) for each 
neutrino trajectory $o- 

Step 8: Find |^ 5) (r + Ar)) := U(H^, r, r + Ar)|^ 0) (r)). 

Step 9: Compute flg^o) = H m&tt {n e {r + Ar)) + A*°(V>< 6 \r + Ar,tf ) for each 
neutrino trajectory i? - 

Step 10: Find \^\r + Ar)) := I7(#bJ> r + ±Ar, iAr)|^ 4) (r + §Ar)). 

Step 11: Find ( r + Ar)) as an average of IV^fy + Ar)) and |'0> 6 ' ) ( r + Ar)) for 
each neutrino mode A. 

The last part of the algorithm estimates the numerical error e and determines the 
course of subsequent computation: 

Step 12: Estimate the numerical error e to be ^ of the maximum of the differences 

(3) (7) 

between \i(j x (r + Ar)) and \ip x (r + Ar)) in all neutrino modes. 

Step 13: Set the new step size to be Ar := £,\/to/£, where £ is an empirical constant, 
and eo is the prescribed error tolerance. 

Step 14: Set \ip\{r + Ar)) := \ipx\ r + Ar )> if 6 < e o- Return to Step 1 and repeat 
the algorithm as necessary. 

Module NBGroup_MM_SMP is a parallelized version of NBGroup_MM using Pthreads. 
In NBGroup_MM_SMP, the algorithm described above is rearranged into a multi-stage 
pipeline through which all neutrino states must flow. Each working thread first enters a 
thread- mutual-exclusive routine NBGroup_PM: :new_task() . In this routine, the thread 
is assigned a subset of neutrino states |V>a) at a particular stage. After finishing 
the operations on at this stage, the thread re-enters NBGroup_PM: :new_task(). 
Subsequently, the thread will be reassigned with another subset consisting of any 
remaining neutrino states at the same stage. If there are no remaining states, and 
if the dependencies are clear, NBGroup_PM: :new_task() will assign to the working 
thread a subset of neutrino states at the next stage in the calculation. The stages are 
arranged so that there is no dependence between neighboring stages except the last 
two. This is done as follows: 

Stage 1: Step 6 in the above algorithm. The last working thread at this stage 
completes step 7. 

Stage 2: Step 2. The last thread also completes step 3. 

Stage 3: Step 8. The last thread also completes step 9. 

Stage 4: Stpdf 4 and 5. 
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Stage 5: Stpdf 10 and 11. 
Stage 6: Stpdf 12, 13, 14 and 1. 

FLAT has been designed to take advantage of computing nodes which are 
capable of running multiple threads simultaneously, e.g., those nodes with multi-cores, 
multiple-CPU's and/or Simultaneous MultiThrcading (SMT) technology. Multiple 
threads can be exploited by FLAT by swapping out the NBGroup_MM module and 
replacing it with NBGroup_MM_SMP. Because Pthreads share a common memory space 
and do not need to pass messages to each other, the parallelization of FLAT using 
NBGroup_MM_SMP is quite efficient. On an IBM 32-way Power4 computing node, multi- 
angle code with 32 working threads can run as fast as ~ 31. 2 x the speed of the 
single-threaded version. On a single-CPU Pentium4 desktop computer with Hyper- 
Threading technology (HT) enabled, threaded single-angle calculations can run > 30% 
faster than non-threaded calculations. 

In order to utilize multiple computing nodes, FLAT is also parallelized using 
Message Passing Interface (MPI) [44] in a manner similar to the implementation in 
BULB (see Section 2.5). This has been done by rewriting portions of the NBGroup_EM 
and NBGroup_I0 modules. For example, function NBGroup_EM: : set_bg() of the 
NBGroup_PNS2D module can be modified to compute the neutrino self-coupling H yu 
based on not only the neutrino states on the residing node but also the states on other 
nodes. Because NBGroup_PM modules are untouched in this scheme, an MPI/Pthread 
hybrid parallel model can be realized by choosing the NBGroup_MM_SMP module together 
with the MPI version of the NBGroup_PNS2D and NBGroup_NETCDF modules. 

2.5. The structure of BULB 

We have developed a second code called BULB to independently verify the results of 
the first code. BULB, developed at LANL, is written in Fortran90 and uses the MPI 
library to communicate between processes. As such it is suitable for use on many 
modern large-scale homogeneous computers. 

BULB is parallelized by splitting neutrino trajectories into different groups of 
emission angles. In this implementation, each node independently handles all neutrino 
energy variables for each angle group. This is acceptable for machines of up to 500- 
1000 nodes, as the largest calculations can be split into one angle per processor. 
For larger machines it would also be desirable to split different energy bins among 
different processors. This may be needed in the future for more complex, non-spherical 
geometries where more than roughly 500 energy bins and 500 angle bins may be 
required. 

The numerical integration algorithm on each node is a fairly simple treatment of 
nonlinear coupled differential equations. For each step, we choose an initial step size 
Ar, evaluating the matter background potential -ff ma tt as an average over the initial 
and final point. In analogy to the procedure described above for FLAT, this increment 
in radius can be translated into the increment in the Affinc parameter At along each 
individual trajectory. 

Initially, we assume that the background neutrino potential is constant over the 
step. The overall Hamiltonian H can then be diagonalized efficiently for 2 x 2 or 
3x3 neutrino cases, and the propagated wave function at the end of the step is 
V'fln = exp(— iiyAi)^>j n i (see Section 2.3). When the neutrino-neutrino scattering 
contribution is small, this formulation can allow for significant propagation distances. 
The step sizes in this case are bounded by the scale height of the matter profile or the 
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length scales related to the fluctuations in the background electron potential. This 
formulation also preserves the unitarity of the neutrino density matrix. 

Upon completion of this initial propagation, the final wave function is used 
to compute a revised background neutrino potential H vv at the final point. This 
new potential is then averaged with the initial potential to produce a neutrino 
potential averaged over the propagation step, similar to the treatment of the electron 
background potential. An improved final wave function is obtained by propagating 
with this averaged interaction. In principle this averaging scheme can be iterated, 
though in practice it is nearly always converged after the first iteration. 

As in FLAT, in BULB we check for convergence with respect to the step size by 
comparing the final density matrix obtained in a single step Sr with that obtained by 
taking two stpdf with half the step size. If the density matrix does not match within a 
predefined tolerance, the step size is reduced and the process is repeated. Periodically, 
we attempt to increase the step size in order to efficiently treat the large radius regions 
where the neutrino background potential and the electron density are fairly small. 

A clear limitation in all these problems is that for each step one must calculate 
the full neutrino density matrix and the angular dependence of the neutrino potential 
for each energy and angle. This involves a sum over all energies, which in the present 
version of BULB is performed independently on each node. Calculating neutrino 
potentials also involves a sum over angles, as the neutrino interaction depends upon 
cos "dpq, where $ pq is the angle between two neutrinos with momenta p and q>, 
respectively. 

For spherical symmetry the angular dependence can be retained by evaluating 
the following sums on each node: 

Pn= £Z>(M) (24a) 

j=l i=l 

j=l i=l 

where index n labels the node, flj is the angle between the j'th neutrino trajectory and 
the radial direction at the current radius, and the sum i and j run over the number of 
energy grid points N e and the number of angles on the local node iVthi respectively. 
The full neutrino potential for all angles can then be evaluated by using MPI global 
sum calls to sum both the angle independent terms p and the angle-dependent terms 
p' and by using the formula 

cos "dij — cos di cos -dj + sin sin dj cos(</>; — (j>j ) , (25) 

where <pi and <fij are the azimuthal angles of the two beams. In spherical symmetry 
the second term vanishes because of the averaging over ((f>i — tfij) [26]. 

Convergence of the overall calculation is carefully checked by comparing results 
with different error tolerances and their associated step sizes, and by comparing results 
with different numbers of energy and angle bins. Although qualitative results can often 
be obtained with a modest number of angular bins, a rather large number of energy 
and angle bins are required to get stable solutions. This can typically involve several 
hundred thousand coupled non-linear differential equations. 

A "typical" calculation such as those described in the results section below involve 
~ (800 angle bins) x (400 energy bins) . Near the surface of the proto-neutron star, the 
wavelengths are very small and integration step sizes can be as small as one millimeter. 
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By the time one reaches 20 km, or about twice the radius of the star, the relevant 
wavelengths have increased and the step size is of order one centimeter. The total 
computational time required for a single calculation is approximately 8000 CPU hours 
on a typical cluster with 2 GHz AMD Opteron processors. The computational time is 
roughly equally divided between computing the evolution on each node and gathering 
the neutrino densities across all nodes. Storage of flavor states for 1000 neutrino modes 
at various radii for later analysis requires approximately 4 GB of storage space. 

For understanding the propagation of the neutrinos from small to large distances, 
BULB can also dump out the neutrino potentials or the full wave function at user- 
defined spatial intervals. Storing the full wave function is rather slow, so this is done 
less often than sums over the wave function or density matrix. The latter arc sufficient 
to determine the rough survival probabilities that could be detected in a terrestrial 
observation. 



3. Results 

The first simulations carried out with both codes focused on the late time, hot bubble 
epoch of the core collapse supernova phenomenon. In this post-explosion environment, 
the matter density surrounding the proto-neutron star is high, but still lower than 
in the earlier shock re-heating epoch. In the hot bubble, the neutrino self-coupling 
Hamiltonian H vv can be dominant in some regions and sub-dominant relative to the 
matter potential in others. 

In our calculations we can adopt a variety of matter density profiles, neutrino 
and antineutrino luminosities, and neutrino and antincutrino energy spectra. For the 
specific example problem of late-time supernova neutrino flavor evolution, we adopt a 
simple analytical profile for electron number density: 

v56„ -3x v _ / Rv-r \ , ,„ n ^ ln 30„ -3s^ /10km x 



n e ~ (1.6 x 10* cm- J )F e cxp 18km j + (6-0 x 10 dU cm- J )V e ^— — J , (26) 

where electron fraction is Y e = 0.4 and the radius of the neutrino sphere is R v = 11 
km. We assume that initial energy spectra for neutrinos are of two-parameter black 
body form 

1 1 El 
MEv) = FM T? cx P (E„/tJ- Vv ) + V (27) 
where r\ v is the degeneracy parameter, T v is the neutrino temperature, and 

x k dx 
exp(x — rj) + 1 

For the late-time supernova environment we take (E Ve ) = 11 McV, (Ep e ) = 16 MeV, 
(E„ T ) = {E 9t ) = 25 MeV, and r] Ve = r) De = tj Vt = r/p T = 3. Here by v T and v T we 
actually mean linear combinations of neutrinos and antineutrinos of the real \i and r 
flavors, respectively. With these choices we have T„ e ~ 2.76 MeV, T Pe ~ 4.01 McV, 
and T Vt = T Vt ~ 6.26 MeV. 

For our example 2x2 mixing treatment presented below, we have taken the 
effective mixing angle to be 9 — 0.1, and mass-squared difference to be Am 2 = 
±3 x 10~ 3 eV 2 . Here the plus sign defines a so-called "normal neutrino mass hierarchy" , 
and the minus sign defines an "inverted neutrino mass hierarchy". There are six 
parameters in the full 3x3 vacuum mixing problem. These parameters consist of 



W = / : (28) 

Jo 
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two mass-squared differences (Am 2 2 and Amf 3 ), three mixing angles (#12, O13 and 
#23) and the CP- violating phase S. Mass-squared differences Am 2 2 and |Amf 3 | and 
mixing angles #12 and 823 have been measured outright, while S and the neutrino mass 
hierarchy are as yet unmeasured. In our 2x2 mixing treatment Am 2 ~ Am 2 3 and 
9 ~ 613. Although 9i3 also remains unmeasured, current laboratory results suggest 
that sin 2 (26»i 3 ) < 0.1. 

Our preliminary single-angle 2-flavor simulations gave very puzzling results. 
Contrary to the studies based on the MSW effect (e.g., [45, 46, 47, 48]), these 
calculations showed that the flavor evolution histories of neutrinos with different 
energies could be coupled together and that all neutrinos and antineutrinos could 
experience collective flavor transformation. This could be true in our calculations 
even when the neutrino-neutrino forward scattering potential was sub-dominant. 

Analysis showed that the collective flavor transformation observed in our 
simulations was not the well known synchronized flavor transformation seen in 
earlier single-angle calculations (e.g., [37, 38]). In synchronized flavor transformation, 
neutrinos and antineutrinos of all energies can undergo a common MSW-like flavor 
transformation in the region where a neutrino with a representative energy would 
experience a conventional MSW resonance. In addition, for neutrino- antineutrino 
gases that are dominated by neutrinos, synchronized flavor transformation can occur 
only in the normal neutrino mass hierarchy case. In our single-angle calculations, 
however, collective flavor transformation occurred for both the normal and inverted 
neutrino mass hierarchies and exhibited oscillatory behaviors over a relatively wide 
radius range (see, e.g., Figure 8 in [26]). 

Based on these preliminary results, we re-investigated collective flavor 
transformation in dense neutrino-antineutrino gases in the presence of ordinary matter 
[24] . We found that the oscillatory collective flavor transformation could be explained 
as "bipolar oscillations" . This type of collective flavor transformation was first studied 
in the context of the early universe [49, 34, 50]. 

Contrary to earlier widely-accepted conclusions, our calculations and analyses 
showed that collective flavor transformation of the bipolar type can appear even in the 
presence of a large and dominant matter density. This is a surprising and paradigm- 
changing result. In the past the astrophysics community had assumed that neutrino 
flavor transformation in the high matter density environment near the neutron star 
would be suppressed (e.g., [37]). However, our results now compelled us to conclude 
that large scale conversion of neutrino flavors could take place in such an environment. 
This conclusion was particularly surprising given the fact that the experimentally- 
inferred neutrino mass-squared differences are small. 

The neutrino and antineutrino trajectories in our calculations are labeled by 
emission angle i?o- This is the angle between the trajectory and the normal to the 
neutrino sphere surface at the neutrino's or antineutrino's point of origin on this 
surface. For example, $0 = corresponds to a radially-directed neutrino trajectory, 
while -do = 7r/2 corresponds to a neutrino moving tangentially to the neutrino sphere. 

In the single-angle approximation, neutrinos and antineutrinos are assumed to 
have the same flavor evolution histories as those with the same energies but traveling 
along radial trajectories. To test the validity of this approximation, we carried out 
multi-angle calculations under the same conditions used in the single-angle cases, 
except of course for the obvious geometric issues associated with binning the trajectory 
directions [26, 51]. 

We first consider the artificial case where we ignore neutrino self-coupling and 
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consider only matter-driven neutrino and antineutrino flavor evolution. This is 
achieved by setting neutrino luminosities L„ c = L j?c = L Vt = L Pt = erg/s. This 
case corresponds to the "MSW" case which had been the paradigm for many if not 
most studies of neutrino flavor transformation effects for the supernova neutrino signal 
(e.g., [45, 46, 47, 48]) and for nucleosynthesis (e.g., [52, 53, 54]). As we will see, this 
case fails to capture the essential physics of the problem. 

The matter density is assumed here to be spherically symmetric. As a result, we 
expect the pure matter-driven MSW case to exhibit little dependence on neutrino or 
antineutrino emission angle (trajectory direction). Figure 3 shows a snapshot (r ~ 100 
km) of one of our simulations for this case. A movie presenting the full simulation 
is available for downloading. There are several prominent features of neutrino and 
antineutrino flavor transformation in this scenario. We see no transformation of 
antineutrinos in this case, and neutrinos are transformed only at relatively low energy 
E v . This is classic MSW behavior: because the neutrino mass-squared difference is 
small, we must go to large radius before neutrinos with average energies are affected, 
and antineutrinos never are. And, indeed, there is little emission angle-dependence in 
this case. 

Figure 4 shows a snapshot (r ~ 100 km) of one of our simulations for the pure 
matter-driven MSW case, but now for the inverted neutrino mass hierarchy. A movie 
presenting the full simulation is available for downloading. These results are in some 
sense similar to those for the normal neutrino mass hierarchy case, except that here 
it is the low energy antineutrino flavors that are transformed. There is no neutrino 
flavor conversion in the inverted neutrino mass hierarchy case. 

We can repeat these calculations for the same supernova model, but now with 
the full neutrino interaction Hamiltonian, including both the matter potential and the 
neutrino self-coupling potential H vv . The resulting neutrino and antineutrino flavor 
evolution is dramatically different from that in the MSW cases. 

Figure 5 shows a snapshot (r ~ 90 km) of one of our simulations for the 
normal neutrino mass hierarchy case but now with the full neutrino self-coupling 
Hamiltonian. A movie presenting the full simulation is available for downloading. In 
this simulation we have taken all neutrino species to have equal luminosities and 
L Ve = Lu e = L Ur = Lp T = 10 51 erg/s. This simulation shows that at values 
of radius r ~ 90 km both neutrinos and antineutrinos over broad energy ranges 
experience large-scale simultaneous flavor transformation. Although earlier analytical 
studies [39, 24] suggested that this behavior was possible, this is nevertheless a 
surprising result as previous numerical simulations did not find this behavior in 
comparable or earlier epochs [37, 38]. In fact, collective flavor transformation in the 
supernova environment previously was thought to be of the synchronized, collective 
form [37]. However, [55] shows that the collective flavor transformation observed in 
this simulation is not synchronization, but rather a neutrino-background-enhanced 
MSW-like flavor transformation. This phenomenon occurs where a neutrino with an 
energy representative of the neutrino-antineutrino gas would encounter a conventional 
MSW resonance. Because the representative energy of a dense neutrino-antineutrino 
gas is usually smaller than the average energy of each neutrino species, this type of 
collective flavor oscillation can occur deeper in the supernova envelope than one would 
expect for pure matter-driven MSW. 

In our simulations we find that neutrinos and antineutrinos emitted in different 
directions behave differently. Of course, this is in direct contradiction to the single- 
angle approximation. In our normal neutrino mass hierarchy simulations, large- 
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Figure 3. For the pure matter-driven MSW case with the normal neutrino 
mass hierarchy we show neutrino survival probability P„„ (upper left panel) and 
antineutrino survival probability Pau (upper right panel) as functions of neutrino 
energy E u or antineutrino energy Ei, and cosine of the emission angle, cosi?o, for 
values of r§ (radius in units of 10 6 cm). The progress bar at the bottom shows 
r§. Survival probability is given by the color code at the top of the panels; red 
corresponding to little or no flavor transformation, and blue corresponding to 
complete flavor conversion. Note that cos$o is unity for radially-propagating 
particles, and zero for tangentially-propagating ones. In the bottom panels, 
the corresponding angle-averaged energy distribution functions / (arbitrary 
normalization) are shown for neutrinos (lower left panel) or antincutrinos (lower 
right panel) for flavor e (t) by the red (blue) curves. The dashed and solid curves 
in the bottom panels are for the neutrino energy spectra at the neutrino sphere and 
at the radius indicated by the progress bar, respectively. This figure corresponds 
to the snapshot at re ~ 10 in the full simulation. 



scale flavor transformation is initiated by neutrinos and antineutrinos in the most 
tangentially-propagating angle bins (i.e., those with low values of cosi?o)- As the 
radius increases, we can see from the full simulation that flavor conversion spreads 
out to other angle bins, eventually causing significant flavor transformation of both 
neutrinos and antineutrinos for nearly all energies and almost all trajectory directions. 

The spreading of significant flavor transformation from tangential to radial 
directions has its origin in the structure of the weak current. More tangentially- 
directed neutrinos or antineutrinos which forward-scatter on the more radially-directed 
background neutrinos and antineutrinos bring larger flavor diagonal and off-diagonal 
potentials to the full self-coupling Hamiltonian H^^. This is a consequence of the 
intersection angle dependence of the forward neutrino-neutrino scattering amplitude 




Figure 4. This figure is the same as Figure 3, except now we show neutrino 
and antineutrino flavor transformation results for the pure matter-driven MSW 
case with the inverted neutrino mass hierarchy. This figure corresponds to the 
snapshot at re ~ 10 in the full simulation. 



[23, 39, 26, 51]. 

As we move out to even larger radius in these simulations, neutrinos and 
antineutrinos experience collective bipolar oscillations. This behavior is evident in the 
full simulation. At large enough radius the neutrino background becomes ineffective 
in influencing flavor transformation. The subsequent neutrino and antineutrino 
flavor evolution is akin to MSW or vacuum neutrino oscillations. However, as will 
be discussed below, the breakdown of collective flavor transformation can leave an 
important "fossil" imprint of the nonlinear neutrino background. 

In the inverted mass hierarchy, neutrino and antineutrino flavor evolution with 
self-coupling is different from the behavior in both the matter-driven MSW case and 
the normal neutrino mass hierarchy self-coupled case. Figure 6 shows a snapshot 
(r ~ 75 km) of one of our simulations for the inverted neutrino mass hierarchy case 
with the full neutrino self-coupling Hamiltonian (L„ e = L Pe = L Vt = L„ T = 10 51 
erg/s). A movie presenting the full simulation is available for downloading. Like the 
normal mass hierarchy case, both neutrinos and antineutrinos can experience large- 
scale, and simultaneous flavor conversion over broad ranges of energies and for nearly 
all emission angles. However, in contrast to the normal mass hierarchy behavior, 
neutrinos and antineutrinos in the inverted mass hierarchy case do not experience 
MSW- like flavor transformation, but rather enter the bipolar oscillation mode directly 
at low radius. Subsequently, as the radius increases, more and more horizontal 




Figure 5. Neutrino and antineutrino flavor transformation results for the 
normal neutrino mass hierarchy case but now with full neutrino self-coupling 
and trajectory coupling. Setups and definitions are the same as in Figure 3. Here 
we take the luminosities of all neutrino species to be 10 51 erg/s. This fi gurc 
corresponds to a snapshot at re, ~ 9.0 in the full simulation. 



fringes appear in the survival probability P vlJ {cosdo : E l y) and Ppp(cos^o,^)- These 
horizontal fringes also appear, although to a less "violent" extent, in the normal mass 
hierarchy case. This phenomenon is associated with multi-angle effects. The one-by- 
one excitation of multipoles in cos t?o so evident in our full simulation is explained by 
the "kinematic decoherence" of bipolar oscillations [56] . 

Neutrino and antineutrino flavor transformation in the inverted neutrino mass 
hierarchy case sets in at a radius rx where the total neutrino flux drops below a critical 
value. This can occur closer to the neutron star than MSW-like transformation in the 
normal neutrino mass hierarchy case. Because rx is solely determined by neutrino 
fluxes and energy spectra in this case [24], neutrinos and antineutrinos can experience 
large-scale flavor transformation even in earlier epochs, where the matter density is 
much higher. Indeed, a simple classical pendulum analogy shows that the inverted 
mass hierarchy gives rise to an inherently unstable neutrino flavor field, like a pencil 
standing on its tip [57, 58]. 

Although the single-angle approximation may not be accurate in general, the 
calculations based on this assumption do capture some if not most of the qualitative 
features present in multi-angle calculations. This has been shown in our simulations 
[26, 51] as well as in the computations carried out by other groups [59, 60]. 

One of the most interesting behaviors seen in both the single-angle and multi- 
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Figure 6. Neutrino and antineutrino flavor transformation results for the 
inverted neutrino mass hierarchy case but now with full neutrino self-coupling 
and trajectory coupling. Setups and definitions are the same as in Figure 3. Here 
we take take the luminosities of all neutrino species to be L v = 10 51 erg/s. This 
plot corresponds to the snapshot at re, ~ 7.5 in the full simulation. 



angle calculations is stepwise spectral swapping (also known as the spectral split) 
[26, 51]. This is where, e.g., v e with energies larger or smaller than some transition 
energy E s swap energy spectra and fluxes with neutrinos of another flavor, say v T . 
The swapping feature is strikingly obvious in, e.g., the later part (larger radius 
part) of the simulation that corresponds to Figures 5. Note the vertical (i.e., for 
all z?o) im e of demarcation between near-zero survival probability for neutrinos with 
E v < E s ~ 10 MeV and significant survival probability for more energetic neutrinos. 
"Swapping" is an apt description of this phenomenon in the sense that there is nearly 
complete conversion of neutrino flavors across a broad range of neutrino energy E v . 
Whether the swapping of spectra occurs for values of E v above or below E s depends 
on the neutrino mass hierarchy. In the normal (inverted) neutrino mass hierarchy 
the swap occurs for neutrinos with energies satisfying E v < E s (E v > E s ). An 
explanation of this phenomenon in the single-angle approximation context has been 
given by [26, 58, 61, 55, 62]. 

In our simulations with the normal neutrino mass hierarchy, the swap energy E s 
decreases as the effective mixing angle 8 (~ 6*13) is decreased. The situation in the 
inverted neutrino mass hierarchy is different. As alluded to above, there is "instability" 
in flavor transformation in the inverted mass scheme. Indeed, our simulations suggest 
that a swap can be obtained in this scheme even for extremely small values of 9 [63]. 
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The direct detection of stepwise swapping in the energy spectra of supernova 
neutrinos could provide a unique way to determine the neutrino mass hierarchy, even 
if #13 is too small to be measured in conventional terrestrial neutrino oscillation 
experiments [63]. This is an important result that could change future supernova 
neutrino detection strategies. The swap energy is relatively small, E s ~ lOMeV, on 
the order of solar neutrino energies. These facts may suggest that supernova neutrino 
detection schemes should focus on low energy, low flux v e 's at a few seconds post-core- 
bounce. 

We have also run 3-neutrino (i.e., 3x3) flavor mixing simulations, where all three 
neutrino flavors are followed simultaneously and sclf-consistently. Given the wealth of 
unexpected phenomena revealed by the nonlinear 2x2 mixing calculations, it would 
be prudent to investigate whether going to a full 3x3 mixing scheme might change 
things yet again. Sometimes the full 3-neutrino mixing problem should be be reducible 
to two separate 2-flavor problems [64, 46, 65, 66]. However, this is not always possible. 

We have run both single-angle and multi-angle tests which employ the full 3-flavor 
mixing scheme for one particular case: the neutronization burst from an O-Ne-Mg 
core-collapse supernova [67, 68]. The neutronization neutrino burst occurs when the 
supernova bounce shock comes through the neutrino sphere at ~ 10 ms post-core- 
bounce. O-Ne-Mg core collapse supernovae have as progenitors stars with masses in 
the range 8 to 12 M . These are very interesting from a neutrino physics standpoint 
because there is little matter overburden for the shock to transit and, as a result, the 
neutrino background can become important very early on even in the normal neutrino 
mass hierarchy case. For this case, in the normal neutrino mass hierarchy, our 3x3 
numerical calculation shows two stepwise spectral swaps for neutrinos in the vacuum 
mass basis, one on the top of the other. This result seems to suggest that, at least 
for this particular case, the full 3-flavor mixing problem is indeed reducible to two 
separate 2-flavor mixing schemes. The result for the inverted mass hierarchy case is 
similar, except that one of the spectral swaps is forbidden because of the conservation 
of a "lepton number" [57, 61]. 

4. Frontiers 

Though there has been dramatic recent progress in understanding the evolution of the 
neutrino flavor field in the supernova environment, there are many outstanding and 
unresolved issues. As recent multi-dimensional hydrodynamics simulations have shown 
(e.g., [10, 12]), the actual supernova matter density could be quite inhomogeneous. 
In addition, it is possible that the neutrino flux and angular distributions could 
deviate significantly from the idealized spherically symmetric models employed in 
our numerical simulations. Nevertheless we believe that our spherically symmetric 
coherent regime calculations at the very least capture the qualitative features of the 
collective flavor oscillation modes discussed above. But a road map for future work 
should include a plan for checking this assertion. 

One of the most significant limitations in current multi-angle simulations is the 
neutrino bulb model employed in these calculations. In the neutrino bulb model, 
neutrinos are assumed to be emitted isotropically from an infinitely thin layer (i.e., 
the neutrino sphere) on the surface of a spherically symmetric black-body (the neutron 
star). In addition, in current neutrino flavor transformation simulations the matter 
profile outside the neutron star is taken to depend only on the radius. 

This model roughly corresponds to the ID supernova simulations, but with 



Simulating nonlinear neutrino flavor evolution 



23 



nowhere near the sophistication of those calculations. For example, it is well known 
that neutrinos with different flavors and with different energies experience different 
opacities and, therefore, decouple from the matter fields at different radii. Even in 
the context of ID supernova models, the initial neutrino energy spectra will show 
deviations from black-body form. Moreover, the neutrino sphere has a finite thickness. 
This will alter the relationship between length parameters and emission angles on 
neutrino trajectories. 

Recent sophisticated supernova simulations suggest that 3D features might play 
significant roles in how massive stars eventually collapse. Needless to say, neutrino 
emission would not be isotropic in these simulations, and the matter profile obviously 
would not be spherically symmetric either. In order for FLAT to treat the 3D 
supernova features in neutrino flavor evolution, a new NBGroup_EM module would be 
necessary. This module would have to model correctly the physics and geometry 
of the problem. BULB would have to be modified in corresponding fashion. Of 
course, actually implementing such modifications in either code would significantly 
enlarge the dimensionality of the problem space and, once achieved, could increase 
the computation time by several orders of magnitude. 

Flavor evolution of neutrinos emitted from the accretion disk around a massive 
black hole is an interesting example of non-sphericity. Even the simplest accretion 
disk model has a more complicated geometry than the neutrino bulb. In principle, 
this problem could be approached in a manner similar to the fD to 3D upgrade, i.e., 
by building a new NBGroup_EM module to treat the special geometric structure of the 
accretion disk. 

Yet another interesting problem is the flavor evolution of the neutrinos and 
antineutrinos in the early universe. This was the environment where the first 
full nonlinear numerical simulations of neutrino flavor evolution were performed 
[33, 69, 70, 49, 34, 50, 71]. These early explorations revealed several interesting 
features of flavor evolution in dense neutrino gases. These features included, for 
example, synchronized and bipolar collective neutrino oscillations. This work has 
had important implications for models of lepton-degenerate cosmologies (see, e.g., 
[72, 73, 74]). Because the early universe and supernovae have, in some sense, similar 
physical environments, the interesting phenomena recently discovered in the supernova 
case might also occur in the early universe, ft could be worthwhile to re-investigate 
neutrino flavor evolution in the early universe too, for example, to see if bipolar 
neutrino oscillations might have consequences for models for lepton number affected 
nucleosynthesis. 

Following neutrino and antincutrino flavor evolution in astrophysical environ- 
ments is a difficult problem. However, astrophysicists are compelled to solve this 
problem because of the potential leverage neutrino flavor transformation has in super- 
nova and early universe physics, and because neutrino flavor mixing is an experimental 
reality. We believe that the solution of the astrophysical neutrino flavor transformation 
problem could usher in a new era where there is a synergistic coupling of terrestrial 
laboratory-based neutrino physics and astrophysical observations. 
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